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ABSTRACT 

Traditionally, the subject of hydromagnetic equilibrium in neutron stars has been 
addressed in the context of standard magnetohydrodynamics, with matter obeying 
a barotropic equation of state. In this paper we take a step towards a more real- 
istic treatment of the problem by considering neutron stars with interior supcrfluid 
components. In this multifluid model stratification associated with a varying matter 
composition (the relative proton to neutron density fraction) enters as a natural in- 
gredient, leading to a non-barotropic system. After formulating the hydromagnetic 
equilibrium of superfluid/superconducting neutron stars as a perturbation problem, 
we focus on the particular case of a three-fluid system consisting of superfluid neutrons 
and normal protons and electrons. We determine the equilibrium structure of dipolar 
magnetic fields with a mixed poloidal-toroidal composition. We find that, with re- 
spect to barotropic models, stratification has the generic effect of leading to equilibria 
with a higher fraction of magnetic energy stored in the toroidal component. However, 
even in models with strong stratification the poloidal and toroidal components are 
comparable, with the former contributing the bulk of the magnetic energy. 
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1 INTRODUCTION 

Neutron stars exhibit a rich phenomenology, observed in a 
variety of channels. The quality of the associated data has 
improved considerably in the last decade, and we are now 
beginning to make detailed inferences about the complex 
physics associated with these systems. A notable recent suc- 
cess concerns the evidence of superfluidity in the compact 
remnant in the Cass iopeia A, the youngest observed neutro n 
star in the Galaxy (|Page et all 1 20 lit IShternin et al.ll201lT ). 
Nevertheless, we are quite far from a truly quantitative 
understanding of these objects. Issues associated with the, 
sometimes extremely strong, magnetic field are particularly 
vexing. The magnetic field is invoked to explain a range of 
observations, yet we do not have particularly reliable the- 
oretical models. Most notably, despite more than 40 years 
of observations, the mechanism that leads to the observed 
radio pulses remains a puzzle. A similar problem concerns 
the observed magnetar giant flares. In particular, we do not 
have a clear picture of the emission mechanism associated 
with the quasiperiodic oscillations seen in the light curve 
of these events. The observed frequencies seem to fit the 
spectrum of magneto-elastic oscillations associated with 



the magnetic field and the elastic crust, but how exactly do 
these mechanical vibrations give rise to the observed X-ray 
variability? Moving inwards, the physics of the neutron star 
core also remains uncertain. The structure of the interior 
magnetic field depends crucially on the composition and the 
state of matter at extreme densities. It is generally expected 
that mature neutron stars are s ufficiently cold to have 
super fl uid/superconducting cores (iBav m. Pethick fc Pines! 
1 19691 : iGlampedakis. Andersson fc Samuelssonl 201ll ). If 
this is the case, then we need to understand how su- 
perconductivity impacts on the long-term magnetic field 
evolution, for example, by expelling the interior field on an 
astrophysically relevant timescale. Alternatively, we need to 
understand how the expected quantised fluxtubes alter the 
magnetic forces and (perhaps) the structure of the global 
field. These issues depend crucially on the matter com- 
position, and the associated critical temperature/density 
for the transition to superfluidity. Current state-of-the-art 
magnetic star modeling is, however, not yet at this level. 
The stark reality is that we do not even understand 
the relevance of varying matter composition. Existing 
models have, almost exclusively, considered barotropic 
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fluid models, that is, models obeying an equation of state 
relation p(p) between the fluid pressure and den si ty (fo r 



Haskell et al. d2008h: Ciolfi et al.l (2009T): 


Lander & Jones 


(2009 


); Ciolfi, Ferrari & Gualtieril (2010 


); Laskv et al. 


(2011 


); ICiolfi et al.l (1201 lh : references to earlier work can 



be found in iMestell ll999j)). Yet, we know that composition 
variations are important for neutron star dynamics, and it 
is generally expected that there may be an impact on the 
ma gnetic field s t ructu re as well. This has been pointed out 
by lReiseneggen (|2009l ) and, to some extent, stratification 
has been incorporated in the numerical models of Braith- 
waite (see f or ex ample iBraithwaite fc Nordlundl (2006); 
iBraithwaitj |2009h ) albeit in terms of a radial entropy 
gradient for an ideal gas equation of state. 

If we want to be able to consider realistic models for 
mature neutron stars, then we need to be able to construct 
non-barotropic magnetic neutron star models. This is the 
aim of the present work. Working in the framework of per- 
turbation theory, which should be adequate for realistic neu- 
tron stars, we will develop models for the global magnetic 
field that takes full account of the multifluid nature of su- 
perfluid neutron stars and the associated variation in the 
internal composition. 

In this paper we focus on hydromagnetic models 
with a neutron superfluid component while 'ignoring' 
proton superconductivity. The magnetic field is assumed 
to have a dipolar structure while for the stellar mat- 
ter we employ a simple polytropic equation of state 
and a phenomenological expression for the stratification 
asso ciated with the proton fraction. A co mpanion pa- 
per (|Lander. Andersson fc Gl"a mpcdakisl l201lh extends the 
modelling to a non-perturbative framework with a generic 
magnetic field geometry, a broader class of non-barotropic 
equations of state and includes a first discussion of the chal- 
lenging issue of hydromagnetic equilibrium in superconduct- 
ing neutron stars. 



2 THE MULTIFLUID MODEL 

Our neutron star model is a multifluid system consisting of 
neutro ns, protons and electrons (e.g. lAndersson fc Corned 
(2006)). The system is assumed to be axisymmetric with the 
fluids at rest. This last assumption is reasonable for slowly 
spinning neutron stars like the magnetars where the mag- 
netic fi eld energy is much larger th an the rotational kinetic 
energy (|Duncan fc Thompson 1992). This work is primarily 
focused on those systems. 

In addition, we make the following assumptions: We 
work in the context of Newtonian gravity which is justi- 
fied given the overall level of precision of our calculation. 
We also assume the presence of neutron superfluidity in 
the entire stellar volume. Real systems, due to the den- 
sity dependence of the critical temperatures for superfluid- 
ity/superconductivity, will have distinct super- and normal 
fluid regions. We ignore the presence of the elastic crust. 
This is an appropriate approximation as long as we focus 
on equilibrium configurations and the crust is in a relaxed 
state. 



2.1 Basic formalism 

The hydromagnetic equilibrium is gov- 
erned by two coupled Euler equat ions (see 
Glampeda kis. Andersson fc Samuelssonl (|201ll) for de- 
tails) : 

V(p x + $) = — Fx, x={n,p} (1) 

px 

In terms of the particle chemical potentials p n ,p P and p, c 
and the mass m — m n — m p we have defined 



fin 



fin 



m 



(2) 



The Euler equations also feature the fluid densities p n ,p P , 
the gravitational potential <E> and the magnetic forces F x (to 
be specified below). The to tal fluid pr essure is given by the 
thermodynamical relation (|Prix||2004 ) 



Vp = PnV/^n + PpVpp 



(3) 



Our approach to the problem of hydromagnetic equilib- 
rium is based on treating the magnetised system as a per- 
turbation about a non-magnetic and spherically symmetric 
"background" configuration. That is, we write /i x — >• p x +<5p; x 
etc. and work at leading order with respect to the per- 
turbations. This approximation implies that the magnetic 
forces should only appear in the perturbed part of equa- 
tion |T} . As an additional simplification we will assume that 
the star retains its spherical shape even in the presence of 
the magnetic field. This is an accurate approximation even 
for the strongest ob served magnetic fields (for example, see 
lHaskell et aTT(l2008l )). 

To begin with, it is easy to show that the non-magnetic 
background configuration is in chemical beta-equilibrium, 
that is, /i n = p, p = p. This follows from eqn. {TJ after setting 
F x = 0, and taking the difference of the two equations. From 
(J3| it follows that 

V8p = p a V6p, n + p P VSp p + SpVp, (4) 

The background pressure is then simply 

Vp = pV/t (5) 

where p — p n + p p is the total density. 

In formulating the perturbation equations it is conve- 
nient to use the perturbed enthalpy 



Sh = — = (1 - x p )8p,„ + x p 8p 
P 

where we have introduced the proton fraction 
„ - Pp 



(6) 



(7) 



We then have 

PnV5p, n + Pp V8p p = p (VSh - SpVxp) (8) 

where 8/3 — 8p, p — Sji n represents the departure from chem- 
ical equilibrium induced by the perturbations. 

The perturbed magnetic equilibrium is then determined 
by the Euler equations 



V(5£x + 5<&) = -F x , x={n,p} (9) 

Px 
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Note that the densities p x appearing here refer to the non- 
magnetic system. These equations are supplemented by the 
Poisson equation for the gravitational potential, 



V 2 <5$ = 4ttGS p 



(10) 



Following a strate gy familiar from the two-fluid os- 
cillati on problem (e.g. lAndersson. Glampedakis fc Haskelll 
(2009)) we can combine the equations ([9]) and produce 
an equivalent pair of "average" and "difference" equations. 
These are 



V(6h + 5$) - 5l3Vx p 



= — F p - — F n 

Pp Pn 



(Fp + F n ) 



(11) 



(12) 



The first equation is essentially the familiar single-fluid Eu- 
ler equation. Notice, however, that apart from the pressure 
and gravitational forces this equation contains an additional 
force term due to Vx p . This extra term represents the effect 
of the composition stratific ation in the hydromagnet i c equi - 
libri um. As discussed by [Reisenegger fc Goldreichl (|l992T ) 
and lPrix fc Rieutordl |2002T l. in the context of neutron star 
dynamics, this term plays the role of an effective buoyancy, 
providing the restoring force for the family of (?-modes in 
oscillating stratified stars. 

This discussion highlights the important fact that a 
stratified system is non-barotropic, that is, the matter is not 
described by a single parameter equation of state of the form 
p = p(p). In the present multifluid system the chemical po- 
tentials depend on both densities, i.e. are functionals of the 
form p x (p n ,p p ). As a result, the equation of state is of the 
biparametric form p = p(p, (3) or, equivale ntly, p = p (p, x p ) 



j|Reisenegger fc Goldrei'chfl992l ; |Pr"ixll2004l ). 

The difference Euler equation (I12|l is 'unique' to multi- 
fluid systems. In the present context it describes the de- 
parture 5/3 from chemical equilibrium driven by the to- 
tal magnetic force. This heterogeneous magneto-chemical 
balance is the main drivin g agent for magnetic ambipolar 
diffusion in neutron stars jG oldrcich & Reisenegger 1992; 
IGlampedakis. Jones fc Sa muelsson 201lh . The fluid motion 
associated with ambipolar diffusion is ignored in the present 
model. This is equivalent to neglecting nuclear reactions 
involving the particles. Inverting this argument, we would 
expect particle reactions (which should be taking place in 
real neutron stars) to induce fluid flow (this is clear from 
the mass continuity equation associated with each fluid). 
This flow would also induce frictional forces appearing in 
the Euler equation (|12|l . The physics of ambipolar flow in 
superfluid neutron stars was recently disc ussed in detail by 
IGlampedakis. Jones fc SamuelssonI ( 201lf ) , who showed that 
the hydromagnetic quasi-equilibrium is still well described 
by {TT} and fl 1 2 p despite the presence of ambipolar diffusion. 
In th e terminology of IGlampedakis. Jones fc Samuelsson] 
l|201ll) this is a 'tranfusion-dominated' quasi-equilibrium. 



2.2 The magnetic forces 

The magnetic forces F x entering in the multifluid formalism 
are directly linked to the properties of neutron star matter, 
specifically to proton superconductivity and neutron super- 
fluidity. 



Neutron stars are expected to manifest both these prop- 
erties in the bulk of their interiors provided their temper- 
ature lies below the critical temperatures T c ~ 10 9 K for 
the onset of proton and neutron superfluidity. According to 
standard cooling theor y this should hap pen roughly a year 
after the star is born (|Page e t al. 2004). Theory also sug- 
gests that the proton superc onductivity is of the type II 
|Bavm. Pethick fc PinesH l969). which means that the mag- 
netic field penetrates the matter by forming a large number 
of quantised magnetic (proton) fluxtubes. 

The magnetohydrodynamics of a multifluid sys- 
tem with a superconducting component is known 
to be significantly different, both qualitatively and 
quant itatively, from s tanda rd magne t ohydro d ynam- 
ics jEasson fc Pethick! Il977l; iMendelj Il99ll . Il998l ; 



IGlampedakis. Andersson fc SamuelssonI 2011 ) . A fun- 



damental difference is the magnetic force itself. 

In the presence of type II proton superconductivity the 
magnetic force exerted on the proton-electron plasma (de- 
noted as F p in the previous Section) is no longer given by the 
familiar Lorentz force. Instead, the superconducting force 
originates from the tension (energy per unit length) asso- 
ciated with the proton fluxtube array. Its ex plicit form is 
l|Glampedakis. Andersson fc Samuelssonll201ll ): 



47T 



B x (V x Hci) + p P V B 



dp P 



(13) 



where B is the smooth-averaged magnetic field and 
Hci = {H cl /B)B. The critical field H cl is directly re- 
lated to the energy per unit length £ p of each individ- 
ual proton fluxtube and is given by (jTillev fc TillevI IT990I ; 
IGlampedakis. Andersson fc SamuelssonI l201ll ) 



4ty 



Pp 



(14) 



where (j>o = hc/2e is the magnetic flux quantum associ- 
ated with a fluxtube. The definition of h c « constant and 
of the entrainment parameter £ <r (p n iPp) can be found in 
IGlampedakis. Andersson fc SamuelssonI (|201ll ). 

Somewhat counter-intuitively, the neutron superfluid 
also experiences a magnetic force F n . This force originates 
from the coupling of the neutron fluid to the proton flux- 
tubes via the entrainment parameter e*, and is given by 
l|Glampedakis. Andersson fc Samu elsson 201lh 



Pn 
4% 



V B 



. dH cl 

dp n 



(15) 



It is worth noting that, with respect to their magnitude, 
F n ~ F p which indicates that both forces are equally impor- 
tant in the equilibrium of superconducting/superfluid neu- 
tron stars. 

A second scenario to be considered is the state where 
there is proton superconductivity without the simultaneous 
presence of neutron superfluidity. This could be an astro- 
physically relevant case if the neutron star core were to be 
relatively hot, with a temperature T somewhere in the range 
5 x 10 s K < T < 10 9 K. This temperature range is identi- 
fied by the recent cooling observations of the neutron star 
located in Ca ssiopeia A as the range for the onset of neutron 
superfluidity l|Page et al.ll201ll ; IShternin et alj|201ll ). In this 
scenario the magnetic force exerted on the proton fluid is 
again given by ()13[) ; this time, however, the entrainment ef- 
fect between neutrons and protons on the mesoscopic scale 
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of individual fluxtubes is absent. This leads to e* = 1 in 
()14[) and the decoupling of the neutron fluid and the mag- 
netic field, i.e. F n = 0. 

A third possible state of matter is that of a sys- 
tem consisting of superfluid neutrons and normal (non- 
superconducting) protons. This could be realised in neutron 
stars provided the bulk magnetic field strength in the stellar 
interior exceeds the critical threshold H C 2 « 10 16 G above 
which superconductivity is suppress ed regardless of temper- 
ature <|Bavm. Pethick fc Pineslll969l l. With superconductiv- 
ity absent, the magnetic field behaves "classically" and F p 
can be identified with the usual Lorentz force Fl, 



F d = F l 



4tt 



(V x B) x B 



(16) 



At the same time, the neutron fluid is oblivious to the pres- 
ence of the magnetic field, i.e. F n = 0. 



3 EQUILIBRIUM WITHOUT 
SUPERCONDUCTIVITY 

Out of the three possible combinations of neutron superflu- 
idity/proton superconductivity the first is likely to be the 
most relevant for neutron stars (and in particular for mag- 
netars). Based on the present understanding of the involved 
physics, we expect the simultaneous presence of superfluid- 
ity and superconductivity. Nevertheless, we will focus on the 
third scenario, where superconductivity is absent. We have a 
good reason for doing so: the calculation of the full supercon- 
ducting/superfluid hydromagnetic equilibrium represents a 
big leap with respect to the existing work on the subject 
which has been limited t o ordinar^l single-flu i d barotropic 



neutro n star models (e.g . Haskell etaLh 12008) ; Ciol fi et al 



20091) : ICiolfi. Ferrari fc Gualtieril l|201of )~ fLander fc Jones 



2009)). By considering a multifluid system with the pro- 



ton superconductivity "switched-off" we can advance our 
understanding of the role of the stratification in the proton 
fraction considerably. 

The hydromagnetic equilibrium for the chosen multi- 
fluid system is described by 



X7(5h + S$) 



- 5pX7x p = -F L 
P 



Pp 



(17) 



(18) 



plus the Poisson equation (|10l) . Eliminating the magnetic 
force between these two equations and integrating, we arrive 
at the Bernoulli-type law, 



Sh + 8$ — XpSf3 — constant 



(19) 



From these equations we can draw some key conclusions. 
Firstly, the magnetic force forbids the establishment of 
chemical equilibrium (excluding the rather special case of a 
force- free field). This point was also discussed at the end of 
Section [2.11 Secondly, we can always write the ratio Fl/p p 
as a perfect gradient (the same is true for Fl/p only in 
the special case of a uniform proton fraction). This second 



1 An exception is the recent work of lAkeiin fc W asscrman l2008f ) 
on toroidal fields in a single-component type II superconducting 
barotropic neutron star model. 



property suggests that, within the present multifluid model, 
the calculation of the hydromagnetic equilibrium should be 
based on the difference Euler equation (|18p rather than the 
total Euler equation Q17|l. The strategy becomes obvious by 
taking the curl of (|18[) 



(20) 



V x <! — Fl } =0 

.Pp 



This is identical to the equation governing the hydromag- 
netic equilibrium in a single-fluid barotropic star after the 
replacement p — > p p . It should be also emphasized that the 
property (|20|l is characteristic of multifluid systems; it is 
not obeyed by stratified systems described by single-fluid 
magnetohydrodynamics. This fact hints at the possibility 
that hydromagnetic equilibrium in a stratified multifluid sys- 
tem may differ fro m that expected in single-fluid systems 
|Reiseneggerll2009l ) . 

Hence, the calculation of the multifluid equilibrium 
should consist of the same steps as the ones taken in the 
more familiar barotropic problem: (i) the B field is de- 
termined by (|20p and the vanishing azimuthal component 
F£ = (a consequence of axisymmetry) for a given back- 
ground density profile (ii) eqns. (|18[1 . (|19[1 and (|10|) can be 
used to find the remaining unknown functions 5/3, 8h, <5$. In 
this work we will concern ourselves with the first step since 
our main interest is the magnetic field configuration in equi- 
librium. The second step would have been essential if we 
wished to calculate the def ormation of the star 's shape due 
to the magnetic field (as in lHaskell et al. I l|200ct) 1. 

The first step in our analysis is to decompose the mag- 
netic field into poloidal and toroidal components, 



B = B F 



Bn 



Bp = VS x V<p, B T = T\7ip (21) 



where S(r,8) and T(r,9) are axisymmetric "stream func- 
tions" with B ■ VS = (we hereafter adopt standard spher- 
ical coordinates r, <p). The explicit form of the magnetic 
field components is 

B" = ^M- B 9 = — B- = -L- ft (22) 
r smS rsmo rsmtf 

It is easy to see that the magnetic field, when written in the 
form (|21|l . is automatically divergence-free. 

As already mentioned, the assumption of axisymmetry 
requires that Ff = 0. This means that we obtain 



Bp • V(zuBt) = 



VS x VT = 



(23) 



which shows that the two stream functions share the same 
level surfaces, i.e. we can write T = T(S) (this is of cou rse 
a classic result, see IChandrasekhar fc Prendergastl (|l956T 0. 

In terms of the stream functions and after some 
straightforward manipulations the Lorentz force becomes 



F L = — 



VS = AVS 



(24) 



where vj 



Atyzu 2 

r sin 6 is the usual cylindrical radius and 

A» = ro 2 V ■ (vj~ 2 V) (25) 



is the so-called Grad-Shafranov operator. We also note that 
eqn. 1|24[) defines the parameter A. 

Specialising to spherical coordinates we can write 

r 



|f • VS + cot0(0- V5)| (26) 
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It is also worth pointing out the relation between this 
operator and Chandrasekha r's "A5 " operator (e.g. 
IChand rasckhar 'fc Prendergastl (| 19561 )). Momentarily 
switching to cylindrical coordinates w, z, ip we can easily 
show that 

A,S = zu 2 (dlP + ^-d^P + dlPj = zu 2 A 5 P (27) 

where P — S/zu 2 . 

Inserting the Lorentz force (|24|) in (|20[) we arrive at 



(28) 



I — I - V,s' = 

,ft>, 



This is solved by (essentially a Grad-Shafranov equation) 

- = F(S) (29) 
Pp 



where F is an arbitrary function, which is in principle 'user- 
specified'. However, the perturbative approach followed here 
(with the magnetic field treated as a perturbation on a 
non-magnetic background and S ~ O(B)) dictates th at the 
unique choice for F(S) is (see also lCiolfi et~ai1 ||2009h ) 



A 

Pp 



F{S) = c + ciS 



(30) 



with co, ci constants (note that co is allowed to be ~ O(S)). 
Thus, the combination of ([24)) . (f26| and (|30]i leads to our 
'final' equation for the hydromagnetic equilibrium 



-AnXppr 2 sin 2 6 (co + ciS) (31) 



So far we have only discussed the equations pertaining 
to the stellar interior. For the exterior (the magnetosphere) 
we use the commonly adopted assumptions of a perfect vac- 
uum and an irrotational field V x B cx = (the index 'ex' 
denotes an exterior quantity). Combined with axisymmetry, 
the irrotationality condition implies a purely poloidal, force- 
free exterior field. From (|21[1 and (|24[) this means that the 
exterior magnetic field is described by 



T cx = 0, 







(32) 



4 DECOMPOSITION IN MULTIPOLES 

The previous equations can be solved by means of an ex- 
pansion in a ngula r spherical harmonics (see, for example, 
ICiolfi et al.l |2009l )). Before doing this, we note that it is 
more convenient to work with the auxiliary function S de- 
fined by 



S = sinedeS 



(33) 



Decomposing in terms of standard spherical harmonics YJ' 
(and fixing m = for our axisymmetric system), 



S=£>(r)l?(0) 



(34) 



we find that A,S takes the particularly simple form 

A»S= -J2 { a i -^+1)^} sin 9deY? (35) 



where a prime denotes a radial derivative. The expansion 
(|34[) also makes direct contact with the multipolar structure 
of the B field itself. For example, we have 

B r (r, 6) = - Y, ^^-a e (r)Y e %6) (36) 



4.1 Poloidal equilibrium 

In this Section we consider a purely poloidal magnetic field 
(T = 0). After some straightforward algebra, we find that 
(|31[) and ([35} lead to the recurrence relation (with £ > 0) 

(£ - l)Q t C t -i -(£ + 2)Q t+1 C e+1 = 



47rr po 



co {Q2S 2 - 2Q 1 8° l } + ci I QeKe-! ae-i 
— Q e+1 N e+1 att+i — (£ — 2>)QtQi-\Qi-2 
+ (£ + A)Q l+1 Q l+2 Qi+3 ae+3 

where 



£(£+!) 



at, 



f 2 



4£ 2 - 1 

2 



(37) 

(38) 
(39) 
(40) 

and 8" is the usual Kronecker-delta. A closer inspection of 
(|37[1 reveals that the even and odd ^-multipoles decouple 
from each other. 

In the exterior space the situation is much simpler, with 
all multipoles decoupling, 



Ke = i(l — Ql+i - Qi+2) + (£+ l)Qe 



m = (£+ 1) (1 - Q\ - Qt-x) + £Q 2 e+1 



! d 2 aT 



£(£ + l)aT = 



C W = 7 (41) 



dr 2 

where di is a constant. 

The above equations are supplemented by boundary 
conditions at the stellar center and surface. Assuming a 
power-law behaviour near r = we find that the desired 
solution is 



J+ 1 







(42) 



At the surface the field can be smoothly matched to the vac- 
uum exterior provided both functions at, a' t are continuous. 
Using (1411) this requirement leads to the surface condition 

a' e = -—a e , r = R (43) 

We now consider the simplest possible configuration, 
namely, a purely dipolar field (ai 7^ 0, ai>2 = 0), in which 
case the only non-trivial equations (|37[) are: 



( 'i -h 2-Kr 2 p p ( 2c + jr c i a i 



2 / ^ 

C\ + 4irr p p ( co + -C\a 



= {I = 0) 
{I = 2) 



(44) 
(45) 



These are mutually consistent provided ci = 0. Thus, the 
dipole poloidal field is described by a single differential equa- 
tion 

« 2 



T ai = -47rc r x p p 



(46) 
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together with 



where the toroidal term contains the angular integral 



S(r,6) = -\ai(r)sin 2 e, A 



(47) 



4.2 Adding a toroidal component 

We now extend the analysis by adding a toroidal magnetic 
field component. Unlike F(S), the function T(S) cannot be 
uniquely specified within the perturbation scheme. 

A general form for T(S) that compl ies with the s pirit o f 
our perturbation scheme is (also used bv lCiolfi et all feootf )) 



T(S) = -( S 



— -1 

So 



e 



(48) 



where £o and So are constant parameters. The former pa- 
rameter represents the overall toroidal field strength while 
the latter designates the poloidal stream function surface 
which acts as the boundary for the toroidal field (this is 
the role of the step-flunction, O, in (|48[) ~). After some simple 
manipulations (|48[1 leads to 



T dT - r 2 ? 



,J£L 

Vol 



+2 S ,h 



(49) 



A common choice in this problem area, and the one 
we will adopt for the rest of this paper, is to choose So as 
the outermost closed stream surface that is tangent to the 
surface of the star. As a result of this choice the T func- 
tion 1)48 | l ensures that the toroidal field is vanishing at the 
surface (as required by the boundary conditions). Moreover, 
this specific choice for T(S) has the advantage of "placing" 
the toroidal component in the vicinity of the poloidal field's 
'neutral line' (i.e. the locus of points where Bp = 0) ; this 
configuration is likely to be dynamically stable as suggested 
by the work of iMarkev fc Tavled l | 1973ft . 

We again assume a purely dipolar field which obviously 
now has a mixed poloidal-toroidal character. This time, how- 
ever, we find that the general recurrence relation for the cte 
coefficients does not terminate at the first equation (as was 
the case for the purely poloidal field). In other words within 
the present framework of a mixed magnetic field we cannot 
have a purely dipolar geometry. 

We sidestep this difficulty by nevertheless assuming 
a dipolar field, truncating the recurrence relation at the 
first equation and setting ci = 0. This simplification 
is not as radical a s it m ay seem. According to the re- 
sults of ICiolfi et all l|2009l ) we would expect the addition 
of higher odd-order multipoles I = 3, 5 etcetera to have 
only a minor impact on the equilibrium structure of the 
1=1 configuration. The veracity of this assertion is fur- 
ther supported by the close agreement between the results 
of this paper and those obtained in the companion paper 
l|Lander. Andersson fc Gla mpedakis 2011) which considers 
a generic magnetic field geometry. 



4.3 Numerical solution 

Following the approach outlined in the preceding paragraph 
we find that (|46[l is replaced by 

at" - -5-oji = -4irc r 2 Xpp + ^Co^ (50) 



T(r) =^ S in*(s-3s|L + ^)<-> 



1 d8 sin 3 8 < 




Aai 


■ 2 

sm 






So 





Aai 



So 



(51) 



To facilitate the numerical integration of ()50)l (or 1)46 ) we 
introduce a new set of dimensionless parameters. First, we 
define x = r/R which also allows us to write 



M ft \ 



(52) 



where M and R are the stellar mass and radius. Expressed 
in terms of the radial field at the magnetic pole we have 



We can then normalise ai as 



ai = cei - 



0) 



2A 



In terms of the new variables (1501) becomes 



d Oil 2 „ a i 2 r I \ ^ a2 rr l \ 



where 



^-^Z, d = capo, Co = (oR 



(53) 
(54) 

(55) 
(56) 



are all dimensionless parameters. The boundary conditions 
at the stellar center and surface accordingly change to 

dx 



a.\(x — > 0) ~ x 



-(1) = -fii(l) (57) 



Our next task is to simplify the toroidal term I(x). 
The angular integrals in (151)1 are non-vanishing within the 
latitudinal interval 9q (r) < 6 < n — 0o(r) with 8o defined as 



|k<Si| sin do = 1, K = B P R /2So = const. 
We can then carry out the angular integrations, 

T( .c) = -2X&! cose <{ i (2 + sin 2 8o) 



(58) 



\\koh \ (3sin 4 6» + 4sin 2 8 + 8) 



5 

\2 /r 6 



+ -r (k<5i) (5 sin b 8 + 6 sin 8 + 8 sin 8 + 16) 



(59) 



Finally, given the previous normalisations, we define the di- 
mensionless magnetic field components 

.r B T 



a i 



B v 
Bp 



1 da\ 

2x dx 



sin 8 



Br, 2x 



sin 8 ( |kq!i I sin 2 — 1)0 (jnai \ sin 2 8 



noting again that B p is the polar magnetic field strength. 
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Figure 1. The proton fraction p rofile x p (p), eqn. 16411, for 7 = 1 
and 7 = 2 is compared to the iKaminker, Haensel fe Yakovlevl 
20011) fit to a typical PAL equation of state 
Prakash. Lattimer &: Ainsworth Il988f ) . The central density 
is taken to be p c = 10 15 gr/cm 3 . This example shows that it is 
reasonable to model the core proton fraction as linear. 



5 RESULTS: MULTIFLUID 

HYDROMAGNETIC EQUILIBRIUM 

In this Section we construct hydromagnetic equilibria by 
solving eqn. (|55p for the radial function ct\. For the non- 
magnetic background we choose a stellar model with an = 1 
polytropic density profile 



M 



\(ttx) 



(63) 



Given the overall precision of our modelling this choice is a 
good approximation to more realistic equations of state. For 
the proton fraction we use 



0.13 



Pc 



(64) 



where p c = kM/AR? 
constant. This expressio n 



the central density and 7 a 
is m otivated by the work of 
iReiseneeeer fc Goldreichl (|l992h who determined x p for a 
Fermi mixture of non-interacting neutrons, protons and elec- 
trons. Comparing the 7 = 1 form of (|64[) against the profile 
of a typical representative of realistic equations of state we 
find good agreement, see Fig. fTJ It should be also pointed 
out that the prescription f|64[) is valid only in the stellar core 
and not in the region of the crust or the surface. We nev- 
ertheless use it for the entire stellar volume; this should be 
a reasonably accurate approximation (see discussion at the 
end of Section [5. 2p . 

For the purpose of this paper we have introduced the 
phenomenological parameter 7 which controls the relative 
proton fraction to density scale-heights, 

x p _ \_p_ 

Vs;p| 7/o' 

and as such it is an effective measure of the stratification. 
Although a value of 7 significantly different from 1 may not 
be such a good approximation to realistic x p profiles (see 
Fig- [Hi it is nevertheless of interest to experiment with this 
parameter to gain intuition about the effect of stratification 
on hydromagnetic equilibrium. 

The numerical integration of (|55[) with the boundary 
conditions is straightforward. For each calculated hy- 
dromagnetic equilibrium we need to specify three parame- 
ters: the proton fraction power-law 7 in (I64|l . the toroidal 
amplitude £o and the ratio k = B P R 2 /2So (without loss of 



L x = 



(65) 



generality we can set k = 1). The remaining constant do is 
fixed during the integration itself. 

A sample of hydromagnetic equilibria is shown in Fig- 
ures [2] and [3] Let us first discuss Fig. [2] which displays the 
radial profile of the normalised magnetic field h(x, 9). In par- 
ticular we show b r along the polar direction 9 = and b e , W 
in the equatorial plane 9 = 7r/2. We have used £o = 0, 10, 20 
(the first is just a poloidal field) and for each case we have 
chosen three different stellar models: 7=1, which is a model 
with 'canonical' composition stratification, 7 = 2, which 
represents a fiducial strongly-stratified model and finally a 
barotropic model, which formally corresponds to x p — 1. 

The equilibria in Fig. [2] are also shown in Fig. [3] in 
the form of contour plots of the normalised stream func- 
tion S/So, projected on an arbitrary meridional plane. The 
two figures have been arranged such that each panel is in 
one to one correspondence. 

The nine displayed equilibria share some common prop- 
erties. They all feature a single neutral line at the equato- 
rial plane (where B 6 (x n ,ir/2) = 0); when a toroidal field 
is present they all resemble a "twisted-torus" configuration 
with the toroidal field roughly occupying the "hole" in the 
vicinity of the neutral line. From this point of view, the 
structure of our mu ltifluid equilibria is similar to the exist- 
ing barotropic ones dCiolfi et alJ2009l : lLander fc Jonej2009l : 
ICiolfi. Ferrari fc Gualtieril l201oh as well as the numerical 
equilib ria of Braithwaite fc Nordlundl (|2006h : iBraithwaitd 
(2009) (the latter ones represent stable hydromagnetic equi- 
libria of a stellar model obeying an ideal gas equation of 
state - see discussion at the end of Section [ED . 



5.1 The role of stratification 

Is composition stratification (non-uniform x p ) an important 
factor in the hydromagnetic equilibrium of multifluid neu- 
tron stars? 

Based on the evidence provided by the poloidal equi- 
libria in Figs. [2]and[3l the answer is clearly yes. A stronger 
stratification (i.e. a larger 7) pushes the location of the neu- 
tral line inwards. The most extreme example shown here 
(the 7 = 2 equilibrium) has its neutral line at x n m 0.6 as 
compared to the x n « 0.8 of the barotropic equilibrium, c.f. 
the left column panels in Figs. [2] and [3] Intuitively, this be- 
haviour makes sense since the force term in (|17[) due to the 
stratification points inwards (the stratification is stable). 

Another interesting effect of stratification regards the 
relative poloidal field amplitude between the stellar surface 
and center. This is again evident in Fig. f2]by looking at each 
column separately (keeping in mind the changing numeri- 
cal scale). A stronger stratification implies a larger central 
poloidal field Bp relative to its surface value (which is the 
same for all equilibria shown here). 

The addition of a toroidal component invariably moves 
the neutral line outwards. As a result of this, an increasing 
toroidal amplitude £0 results in a smaller portion of the stel- 
lar volume being occupied by the toroidal field. At the same 
time, the overall magnetic configuration appears to take 
similar forms regardless of stratification. For instance, all 
£0 = 20 equi libria look s i milar to the single-fluid barotropic 
equilibria of ICiolfi et ail (|2009h . 

Another way of looking at the impact of stratification 
on the magnetic field structure is by calculating the mag- 
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netic energy stored in the poloidal and toroidal components. 
These are given by 



E P 



1 

8^ 



dVBy 



1 

8^ 



dVB T 



(66) 



where the integrals are taken over the stellar volume. 

A useful quantity to display is the ratio Et/Ep; which 
can easily be expressed in terms of the normalised parame- 
ters of Section 14.21 Our results are shown in Fig. 2] where 
we have again considered the same three models as before 
(barotropic, 7 = 1, 7 = 2). 

The message from Fig. U is rather clear: stratification 
leads to hydromagnetic equilibria with a higher Et /Ep en- 
ergy ratio content with respect to their barotropic counter- 
parts. However, even for our most strongly stratified model 
the actual energy ratio "saturates" at a small value, never 
exceeding the level of Et/Ep ~ 0.1. This is a direct conse- 
quence of the outward "motion" of the neutral line with an 
increasing toroidal component. 

In terms of the actual strength of the toroidal field, we 
note (Fig. [2]) that the maximum relative magnitude \B V /B p \ 
(located at 8 = 7r/2 ) never exceeds a factor ^3 — 5. 

At this point it would be useful to confront the re- 
sults in Figs. [2] and [5] wit h those by Bra i thwaite and collab - 
orators (|Braithwaite fc Nordlundl I2OO6I ; fBraithwaitel [2009). 
This comparison is of particular interest given that Braith- 
waite's stellar model is the only (previous) non-barotropic 
model in the literature. However, this is where the similar- 
ities between the two models end. Braithwaite utilises an 
ideal gas equation of state to build a stellar model charac- 
terised by a stratification in the temperature/entropy pro- 
file. This may be an appropriate model for a main sequence 
star but it is not so good for neutron stars which are highly 
isothermal/isentropic objects with stratification associated 
with matter composition. Despite this fundamental differ- 
ence, the stable quasi-axisymmetric twisted-torus equilibria 
obtained by Braithwaite look roughly the same as the ones 
shown in our Fig. [3] However, Braithwaite's equilibria have 
a much higher energy ratio Et/Ep 3> 1; in that sense they 
are dominated by the toroidal field. In our opinion, the most 
natural explanation for this variance in the magnetic ener- 
gies Et/Ep appears to be the difference in the physics of 
stratification between the two models. An alternative (and 
less likely) explanation could be that the equilibria in Fig. 
are for some reason dynamically unstable and there is a sec- 
ond, yet undiscovered, class of axisymmetric hydromagnetic 
solutions with much stronger toroidal fields. 



5.2 The effect of a surface current 

A different aspect of hydromagnetic equilibrium in neutron 
stars that has received little attention so far is the possibil- 
ity of having electric currents located at the stellar surface. 
Realistic neutron stars could in fact harbour such currents 
as a result of the presence of the magnetosphere. 

A detailed modelling of surface currents in neutron stars 
is beyond the scope of this paper. Indeed, we are agnostic 
on the question of how such currents could be generated 
and maintained. We may nevertheless learn something useful 
by following a simple phenomenological approach where the 
presence of a surface current is mimicked by a modification 
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0.2 0.4 0. 



x = r/R 



Figure 2. The radial profile of the normalised magnetic field 
components b l = B l /B p (B p is the field at the magnetic pole and 
i = {r,8,ip}) for three stellar models: (i) a barotrope (top row), 
(ii) a stratified model with a 7 = 1 profile for the proton fraction 
x p (middle row) and (iii) a stratified model with 7 = 2 (bot- 
tom row). Each column corresponds to a different choice for the 
toroidal field amplitude, £0 = 0, 10, 20. The V field is calculated 
along the 9 = direction while the b B and W components are 
calculated along 8 = n/2. Note that the numerical y-axis scale is 
kept fixed between the columns of a given row. 



of the magnetic field at the stellar surface. A generic con- 
sequence of a surface current is to make the magnetic field 
discontinuous at the surface (this is a direct consequence of 
Ampere's law V x B = 47rJ/c). The idea then is to adopt 
a slightly different set of surface boundary conditions which 
incorporate such a discontinuity. A possible choice i^f] 



h r = bl 



W = (x = l) (67) 



which corresponds to an azimuthal surface current. The 
parameter £ controls the surface current 'strength' (with 
£ = 1 representing the previous case of a continuous B 
field and vanishing surface current). In terms of the ol\ func- 
tion the new set of surface conditions (|67[) is equivalent to 
ddi/dx = — £qi. 

For the purpose of illustration we have constructed two 
equilibria using the stratified 7=1 model with £0 = 10 and 
for two fiducial values £ = 2 and £ = 4 for the discontinu- 
ity parameter. The results are shown in Fig. [5] The most 
notable difference with respect to the previous equilibria re- 
lates to the larger toroidal field component and energy ratio 
Et/Ep; the latter can easily exceed the few percent level. 
Indeed it is now possible to build equilibria with comparable 
energies, i.e. Et ~ Ep (e.g. the £ = 4 model in Fig. GsJ. In a 
sense, the equilibria discussed in this section can be viewed 
as the middle ground between magnetic configurations that 
are f ully confined within the star and can have Et ^> Ep 
(e.g. lHaskell et ah! ((2008)) and the configurations discussed 
in the previous sections which smoothly extend to the ex- 
terior space and have Et <C -Ep. If nothing else, this result 

2 The radial component B r is always continuous at an interface 
as a result of V ■ B = 0. 
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Figure 3. The hydro-magnetic equilibria of Fig. [2] shown as con- 
tour plots of the normalised stream function S/So. Each panel 
corresponds to a panel of Fig. [2] located in the same row and col- 
umn. Top row: barotropic models. Middle row: stratified 7 = 1 
models. Bottom row: stratified 7 = 2 models. Columns from left 
to right: £o = 0, 10, 20. The thick circular line represents the 
stellar surface. The toroidal component b v occupies the region 
bounded by the outermost closed poloidal streamline and it is 
shown with dashed lines. 




Figure 4. The ratio E^/Ep of toroidal to poloidal magnetic 
field energy as a function of the dimensionless toroidal amplitude 
fo- The stellar models are the same as in the previous Figures, 
namely, a barotrope, and 7 = 1, 7 = 2 proton fraction strati- 
fication. For the models shown here the maximum toroidal field 
strength is £?™ ax < 4B P , where B p is the polar magnetic field. 
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6 6 




" 6 0.2 0.4 0.6 0.8 1 0.2 0.4 (1.6 0.8 l" 

x = r/R x = r/R 



Figure 5. Hydromagnetic equilibria for a 7 = 1, £0 = 10 model 
with a "surface current" (see Section [5.2l for details). Each contour 
plot corresponds to the radial profiles plot below it (see Figs. [2] 
and[3]for a description). Left panels: £ = 2. Right panels: £ = 4. 
The toroidal to poloidal energy ratio for these two equilibria is 
E T /E P = 0.267 (left), E T /E P = 0.654 (right). 

highlights the need for a detailed study of the effect that a 
non-vacuum magnetosphere could have on the structure of 
the hydromagnetic equilibrium in realistic neutron stars. 

Before moving on, we note that the truncated recur- 
rence relation also allow a second class of solutions, corre- 
sponding to a different value of the constant do. This new 
family of solutions appears only for mixed poloidal-toroidal 
fields, in both barotropic and stratified models, but as it is 
likely to be an artifact of the truncation in the multipole 
expansion of the stre am f unction S (c . f. the discussion of 
IColaiuda et all (|2008T ) and ICiolfi et al.l (|2009l )) we will not 
consider it in detail here. 

We also carried out a set of computations in order to 
assess the effect of adding a fiducial spherical layer of mat- 
ter with a barotropic (i.e. x p = 1) equation of state. This 
construction can be considered as a toy model of a neutron 
star with a crust. In our model the 'crust' was added to 
the outer x > 0.9 region of the star. Not surprisingly, the 
resulting equilibria only display minor differences from our 
previous results. 

Finally, we carried out computations of equilibria with 
a non-zero constant ci (see eqns. (|31[) and (|44|l 'l. For a broad 
range of ci values (including the case where Co = and ci 
is fixed during the integration itself) the resulting solutions 
display only minor differences with respect to the ones shown 
in Figs. [2] and H 



6 CONCLUDING DISCUSSION 

This paper, together with its companion 
l|Lander. Andersson fc Glampedakisl l201ll) , represents 
the first detailed study of the hydromagnetic equilib- 
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rium in multifluid neutron stars. In spite of the obvious 
simplicity of our models (axisymmetry, dipolar magnetic 
field geometry, omission of proton superconductivity) we 
have managed to probe important aspects of multifluid 
magnetohydrodynamics . 

We have shown that the problem of hydromagnetic 
equilibrium in a multifluid system, consisting of superfluid 
neutrons and normal protons, can be formulated in terms 
of a Grad-Shafranov equation despite the inherent non- 
barotropicity of the system due to the varying composition. 

The main result of this work, summarised in Figs. OH 
concerns the impact of the composition stratification (ex- 
pressed in terms of the proton fraction gradient V:r p ) on 
the structure of the magnetic field in equilibrium. The ob- 
tained equilibria are of the twisted-torus type, containing a 
mixture of poloidal and toroidal magnetic field components. 
We find that stratification generally increases the energy 
stored in the toroidal field by allowing a larger portion of 
the stellar volume to be occupied by it. However, even for 
unrealistically strong stratification (our 7 = 2 model), we 
conclude that the produced equilibria are still dominated by 
the poloidal field, with the toroidal magnetic energy never 
exceeding a fraction ~ 0.1 of the poloidal energy (Fig. 2}. At 
first glance, our results see m at odds with the non-barotropi c 
twisted-toru s equil ibria of iBraithwaite fc Nordlundl (2006); 
iBraithwaitei (|2009h . in the sense that the latter are domi- 
nated by the toroidal field. As discussed in Section 15.11 the 
discrepancy may have a perfectly natural explanation in the 
different type of stratification assumed in the underlying 
stellar models. Nonetheless, this issue needs to be clarified 
in future work. 

The main conclusions of this paper are broadly 
supported by the more de tai led modelling of 
lLander. Andersson fc Glampedakis] (|201lT ) . which ex- 
tends the present work to general, non-dipolar magnetic 
configurations using a wider class of non-barotropic equa- 
tions of state. The same work also lays the basis for 
the modelling of more realistic equilibria with type II 
superconductivity. 

A crucial issue not addressed here is that of the sta- 
bility of the calculated equilibria. Recent numerical sim- 
ulatio ns l|Braithwaitei 120081 : lLaskv et all l201ll ; ICiolfi et al] 
20 1 1 ) have revealed the existence of dynamically stable 
non-axisymmetric equilibria in both barotropic and non- 
barotropic (ideal gas) stellar models. Thus, it is of obvious 
importance to assess the stability of axisymmetric equilibria 
like the ones obtained here. 

Our multifluid hydromagnetic equilibria provide an 
improved understanding of different aspects of magnetar 
physics. Taken at face value our results suggest that the 
surface magnetic field in ma gnetars (e.g. as inferred by the 
standard spin-down formula (|Shapiro fc Teu kolskv 198j|)) is 
also a good indicator of the interior magnetic field strength. 
More specifically, and contrary to a widespread belief in the 
literature, the maximum toroidal field in the stellar inte- 
rior is likely to be comparable to the value at the magnetic 
pole. This result would suggest that the reservoir of mag- 
netic energy available for powering giant flares in magnetars 
is constrained to a level indicated by the surface field. This 
may impact on the viability of some proposed scenarios. 

Finally, the moderate structural differences in the mag- 
netic field geometry caused by the stratification in a mul- 



tifluid system is likely to have some impact on the ongo- 
ing theoretical effort to interpret the quasi-periodic oscil- 
lations obs erved during giant f l ares a s magnetar pulsatio n 
modes (e.g.lvan Hoven fc Levinl l|201 j ): lGabler et al.l(201ll ); 
IColaiuda fc Kokkotasl 1 201ll )V It should be clear that any 
model aiming to confront the real data must incorporate 
the key aspects of multifluid magnetohydrodynamics. 
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